/**
 * 
 */
package ua.lviv.franko.mlsFunctions.twoDimension;

import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

import org.jzy3d.maths.Range;
import org.jzy3d.plot3d.builder.Mapper;

import ua.lviv.franko.integration.Gauss2D;
import ua.lviv.franko.integration.IFunction2D;
import ua.lviv.franko.matrix.SimpleLU;
import ua.lviv.franko.matrix.hash.gauss.Gauss;
import ua.lviv.franko.mlsFunctions.oneDimension.weightFunctions.CubicSpline;
import ua.lviv.franko.mlsFunctions.twoDimension.weightFunctions.IWeightFunction2D;
import ua.lviv.franko.mlsFunctions.twoDimension.weightFunctions.SquareSupport;
import ua.lviv.franko.visualisation.Visualise3D;

//import visualisation.EFGResult2D;

/**
 * @author marko
 * 
 */
public class BaseFunction {
	protected static final Map<Particle, List<Double>>	storage	= new HashMap<Particle, List<Double>>();
	protected int										PolinomialDegree;
	protected Particle									xi;
	protected IWeightFunction2D							func;
	protected double[][]								A;
	protected double[]									b;
	protected List<Particle>							NodeList;
	protected double									EPS		= 0.0000001;

	public BaseFunction(int PolinomialDegree, Particle xi, IWeightFunction2D func) {
		this.PolinomialDegree = PolinomialDegree;
		this.xi = xi;
		this.func = func;
		storage.put(xi, getBaseVector(xi.xi, xi.yi, xi));
	}

	public double calculate(double x, double y) {
		List<Double> baseVector = storage.get(xi);
		if (baseVector == null) {
			baseVector = getBaseVector(xi.xi, xi.yi, xi);
			storage.put(xi, baseVector);
		}
		int n = baseVector.size();
		b = new double[n];
		A = new double[n][n];
		for (int i = 0; i < NodeList.size(); i++) {
			Particle particle = NodeList.get(i);
			calculateAlfa(x, y, particle);
		}

		double val = func.calculate(x, y, xi.xi, xi.yi, xi.rho);
		boolean key = false;
		for (int i = 0; i < b.length; i++) {
			b[i] = baseVector.get(i) * val;
			if (Math.abs(b[i]) > EPS) {
				key = true;
			}
		}

		double returnValue = 0.0;
		if (key) {
			double[] res = Gauss.calculate(A, b);

			baseVector = getBaseVector(x, y, x, y, 0.0);
			for (int i = 0; i < b.length; i++) {
				returnValue += (baseVector.get(i)) * res[i];// returnValue+=BaseVector.get(i)*res[i];
			}
		}
		if (Double.isNaN(returnValue)) {
			return 0;
		}
		return returnValue;
	}

	public double[] calculateAll(double x, double y) {
		if ((xi.xi - xi.rho > x) || (xi.xi + xi.rho < x)) {
			return new double[] { 0, 0, 0 };
		}
		if ((xi.yi - xi.rho > y) || (xi.yi + xi.rho < y)) {
			return new double[] { 0, 0, 0 };
		}
		List<Double> baseVector = storage.get(xi);
		if (baseVector == null) {
			baseVector = getBaseVector(xi.xi, xi.yi, xi);
			storage.put(xi, baseVector);
		}
		int n = baseVector.size();
		b = new double[n];
		// ???????? ??????? ?????? ?
		A = new double[n][n];
		for (int i = 0; i < NodeList.size(); i++) {
			Particle particle = NodeList.get(i);
			calculateAlfa(x, y, particle);
		}

		double val = func.calculate(x, y, xi.xi, xi.yi, xi.rho);
		boolean key = false;
		for (int i = 0; i < b.length; i++) {
			b[i] = baseVector.get(i) * val;
			if (Math.abs(b[i]) > EPS) {
				key = true;
			}
		}
		List<Double> baseVectorSimple = getBaseVector(x, y, x, y, 0.0);
		List<Double> baseVectorSimpleX = getBaseVectorDerivatyX(x, y, x, y, 0.0);
		List<Double> baseVectorSimpleY = getBaseVectorDerivatyY(x, y, x, y, 0.0);
		double returnValue = 0.0;
		double[] ret = new double[3];
		if (key) {
			SimpleLU lu = new SimpleLU(A);
			double[] res = lu.calculate(b);

			for (int i = 0; i < b.length; i++) {
				returnValue += baseVectorSimple.get(i) * res[i];
			}
			if (Double.isNaN(returnValue)) {
				returnValue = 0;
			}
			ret[0] = returnValue;

			// calculate derivaty X
			returnValue = 0.0;
			key = false;
			double[] tmp = new double[b.length];
			double[] tmpY = new double[b.length];
			for (int j = 0; j < this.NodeList.size(); j++) { // ?? ??????
				Particle xj = NodeList.get(j);
				val = func.calculateDerivatyX(x, y, xj.xi, xj.yi, xj.rho);// func.evaluateDerivaty((x-xj.x)/xj.rho);
				double valY = func.calculateDerivatyY(x, y, xj.xi, xj.yi, xj.rho);
				baseVector = storage.get(xj);
				if (baseVector == null) {
					baseVector = getBaseVector(xj.xi, xj.yi, xj);
					storage.put(xj, baseVector);
				}
				double tmpVal = 0.0;
				for (int k = 0; k < res.length; k++) {// ????? ?? ?????? ? ? ???? ...
					tmpVal += baseVector.get(k) * res[k];
				}

				for (int k = 0; k < res.length; k++) {// ????? ?? ?????? ?
					tmp[k] += baseVector.get(k) * tmpVal * val;
					tmpY[k] += baseVector.get(k) * tmpVal * valY;
				}
			}
			val = func.calculateDerivatyX(x, y, xi.xi, xi.yi, xi.rho);// (x, xi.xi, xi.rho); func.evaluateDerivaty((x-xii)/rho);
			baseVector = storage.get(xi);
			for (int i = 0; i < b.length; i++) {
				b[i] = baseVector.get(i) * val - tmp[i];
				if (Math.abs(b[i]) > EPS) {
					key = true;
				}
			}
			for (int j = 0; j < res.length; j++) {
				returnValue += baseVectorSimpleX.get(j) * res[j];
			}
			if (key) {
				double[] res2 = lu.calculate(b);

				for (int j = 0; j < res2.length; j++) {
					returnValue += baseVectorSimple.get(j) * res2[j];
				}
			}
			ret[1] = returnValue;
			// calcelate derivaty y

			// calculate derivaty X
			returnValue = 0.0;
			key = false;
			// tmp = new double[b.length];
			// for (int j = 0; j < this.NodeList.size(); j++) { // ?? ??????
			// Particle xj = NodeList.get(j);
			// val = func.calculateDerivatyY(x,y,xj.xi,xj.yi,xj.rho);//
			// func.evaluateDerivaty((x-xj.x)/xj.rho);
			// BaseVector = this.getBaseVector(xj.xi,xj.yi, xj);
			// double tmpVal = 0.0;
			// for (int k = 0; k < res.length; k++) {// ????? ?? ?????? ? ? ????
			// ...
			// tmpVal += BaseVector.get(k)*res[k];
			// }
			//
			// for (int k = 0; k < res.length; k++) {//????? ?? ?????? ?
			// tmp[k]+=BaseVector.get(k)*tmpVal*val;
			// }
			// }
			val = func.calculateDerivatyY(x, y, xi.xi, xi.yi, xi.rho);
			baseVector = storage.get(xi);
			for (int i = 0; i < b.length; i++) {
				b[i] = baseVector.get(i) * val - tmpY[i];
				if (Math.abs(b[i]) > EPS) {
					key = true;
				}
			}
			for (int j = 0; j < res.length; j++) {
				returnValue += baseVectorSimpleY.get(j) * res[j];
			}
			if (key) {
				double[] res2 = lu.calculate(b);

				for (int j = 0; j < res2.length; j++) {
					returnValue += baseVectorSimple.get(j) * res2[j];
				}
			}
			ret[2] = returnValue;

		}

		return ret;
	}

	public Values[] calculateAllFunctions(double x, double y) {
		List<Double> baseVector = storage.get(xi);
		int n = baseVector.size();
		double[][] bb = new double[n][NodeList.size()];
		// ???????? ??????? ?????? ?
		A = new double[n][n];
		for (int i = 0; i < NodeList.size(); i++) {
			Particle particle = NodeList.get(i);
			calculateAlfa(x, y, particle);
		}

		double val = func.calculate(x, y, xi.xi, xi.yi, xi.rho);
		// BaseVector = getBaseVector(x, new Particle(x, 0.0));
		// boolean key=true;
		for (int i = 0; i < bb[0].length; i++) {
			Particle p = NodeList.get(i);

			baseVector = storage.get(p);
			if (baseVector == null) {
				baseVector = getBaseVector(p.xi, p.yi, p);
				storage.put(p, baseVector);
			}
			val = func.calculate(x, y, p.xi, p.yi, p.rho);
			for (int j = 0; j < bb.length; j++) {
				bb[j][i] = baseVector.get(j) * val;
			}
		}

		Values[] returnValue = new Values[NodeList.size()];
		SimpleLU lu = new SimpleLU(A);
		double[][] res = lu.calculate(bb);

		List<Double> baseVectorSimple = getBaseVector(x, y, x, y, 0.0);
		List<Double> baseVectorSimpleX = getBaseVectorDerivatyX(x, y, x, y, 0.0);
		List<Double> baseVectorSimpleY = getBaseVectorDerivatyY(x, y, x, y, 0.0);
		for (int i = 0; i < bb[0].length; i++) {
			double valF = 0.0;
			double valFx = 0.0;
			double valFy = 0.0;
			Particle p = NodeList.get(i);
			if ((p.xi - p.rho > x) || (p.xi + p.rho < x)) {
				returnValue[i] = new Values(valF, valFx, valFy);
				continue;
			}
			if ((p.yi - p.rho > y) || (p.yi + p.rho < y)) {
				returnValue[i] = new Values(valF, valFx, valFy);
				continue;
			}
			for (int j = 0; j < bb.length; j++) {
				valF += baseVectorSimple.get(j) * res[j][i];
				valFx += baseVectorSimpleX.get(j) * res[j][i];
				valFy += baseVectorSimpleY.get(j) * res[j][i];
			}
			returnValue[i] = new Values(valF, valFx, valFy);
		}

		return returnValue;
	}

	public double calculateDerivatyX(double x, double y) {

		if ((xi.xi - xi.rho > x) || (xi.xi + xi.rho < x)) {
			return 0;
		}
		if ((xi.yi - xi.rho > y) || (xi.yi + xi.rho < y)) {
			return 0;
		}

		ArrayList<Double> BaseVector = getBaseVector(xi.xi, xi.yi, xi);
		int n = BaseVector.size();
		b = new double[n];
		// ???????? ??????? ?????? ?
		A = new double[n][n];
		for (int i = 0; i < NodeList.size(); i++) {
			Particle particle = NodeList.get(i);
			calculateAlfa(x, y, particle);
		}

		double val = func.calculate(x, y, xi.xi, xi.yi, xi.rho);// (x, xi.xi, xi.rho);
		// BaseVector = getBaseVector(x, new Particle(x, 0.0));
		boolean key = false;
		for (int i = 0; i < b.length; i++) {
			b[i] = BaseVector.get(i) * val;
			if (Math.abs(b[i]) > EPS) {
				key = true;
			}
		}

		double returnValue = 0.0;
		if (key) {
			SimpleLU lu = new SimpleLU(A);
			double[] res = lu.calculate(b);

			// val = func.evaluateDerivaty((x-xii)/rho);
			key = false;
			double[] tmp = new double[b.length];
			for (int j = 0; j < this.NodeList.size(); j++) { // ?? ??????
				Particle xj = NodeList.get(j);
				val = func.calculateDerivatyX(x, y, xj.xi, xj.yi, xj.rho);// func.evaluateDerivaty((x-xj.x)/xj.rho);
				BaseVector = this.getBaseVector(xj.xi, xj.yi, xj);
				double tmpVal = 0.0;
				for (int k = 0; k < res.length; k++) {// ????? ?? ?????? ? ?
														// ???? ...
					tmpVal += BaseVector.get(k) * res[k];
				}

				for (int k = 0; k < res.length; k++) {// ????? ?? ?????? ?
					tmp[k] += BaseVector.get(k) * tmpVal * val;
				}
			}
			val = func.calculateDerivatyX(x, y, xi.xi, xi.yi, xi.rho);// (x,
																		// xi.xi,
																		// xi.rho);
																		// //func.evaluateDerivaty((x-xii)/rho);
			BaseVector = this.getBaseVector(xi.xi, xi.yi, xi);
			for (int i = 0; i < b.length; i++) {
				b[i] = BaseVector.get(i) * val - tmp[i];
				if (Math.abs(b[i]) > EPS) {
					key = true;
				}
			}
			BaseVector = this.getBaseVectorDerivatyX(x, y, new Particle(x, y, 0.0));
			for (int j = 0; j < res.length; j++) {
				returnValue += BaseVector.get(j) * res[j];
			}
			if (key) {
				double[] res2 = lu.calculate(b);

				BaseVector = this.getBaseVector(x, y, new Particle(x, y, 0.0));
				for (int j = 0; j < res2.length; j++) {
					returnValue += BaseVector.get(j) * res2[j];
				}
			}

		}
		if (Double.isNaN(returnValue)) {
			return 0;
		}
		return returnValue;
		// double h=0.0000001;
		// return((this.calculate(x, y) - this.calculate(x-h, y) )/h);
		// return 0.0;

	}

	public double calculateDerivatyY(double x, double y) {

		if ((xi.xi - xi.rho > x) || (xi.xi + xi.rho < x)) {
			return 0;
		}
		if ((xi.yi - xi.rho > y) || (xi.yi + xi.rho < y)) {
			return 0;
		}

		ArrayList<Double> BaseVector = getBaseVector(xi.xi, xi.yi, xi);
		int n = BaseVector.size();
		b = new double[n];
		// ???????? ??????? ?????? ?
		A = new double[n][n];
		for (int i = 0; i < NodeList.size(); i++) {
			Particle particle = NodeList.get(i);
			calculateAlfa(x, y, particle);
		}

		double val = func.calculate(x, y, xi.xi, xi.yi, xi.rho);// (x, xi.xi,
																// xi.rho);
		// BaseVector = getBaseVector(x, new Particle(x, 0.0));
		boolean key = false;
		for (int i = 0; i < b.length; i++) {
			b[i] = BaseVector.get(i) * val;
			if (Math.abs(b[i]) > EPS) {
				key = true;
			}
		}

		double returnValue = 0.0;
		if (key) {
			SimpleLU lu = new SimpleLU(A);
			double[] res = lu.calculate(b);

			// val = func.evaluateDerivaty((x-xii)/rho);
			key = false;
			double[] tmp = new double[b.length];
			for (int j = 0; j < this.NodeList.size(); j++) { // ?? ??????
				Particle xj = NodeList.get(j);
				val = func.calculateDerivatyY(x, y, xj.xi, xj.yi, xj.rho);// func.evaluateDerivaty((x-xj.x)/xj.rho);
				BaseVector = this.getBaseVector(xj.xi, xj.yi, xj);
				double tmpVal = 0.0;
				for (int k = 0; k < res.length; k++) {// ????? ?? ?????? ? ?
														// ???? ...
					tmpVal += BaseVector.get(k) * res[k];
				}

				for (int k = 0; k < res.length; k++) {// ????? ?? ?????? ?
					tmp[k] += BaseVector.get(k) * tmpVal * val;
				}
			}
			val = func.calculateDerivatyY(x, y, xi.xi, xi.yi, xi.rho);// (x,
																		// xi.xi,
																		// xi.rho);
																		// //func.evaluateDerivaty((x-xii)/rho);
			BaseVector = this.getBaseVector(xi.xi, xi.yi, xi);
			for (int i = 0; i < b.length; i++) {
				b[i] = BaseVector.get(i) * val - tmp[i];
				if (Math.abs(b[i]) > EPS) {
					key = true;
				}
			}
			BaseVector = this.getBaseVectorDerivatyY(x, y, new Particle(x, y, 0.0));
			for (int j = 0; j < res.length; j++) {
				returnValue += BaseVector.get(j) * res[j];
			}
			if (key) {
				double[] res2 = lu.calculate(b);

				BaseVector = this.getBaseVector(x, y, x, y, 0.0);
				for (int j = 0; j < res2.length; j++) {
					returnValue += BaseVector.get(j) * res2[j];
				}
			}

		}
		if (Double.isNaN(returnValue)) {
			return 0;
		}
		return returnValue;
		// double h = 0.0000001;
		// return((this.calculate(x, y) - this.calculate(x, y-h) )/h);
		// return 0.0;

	}

	public void display() {
		final BaseFunction func = this;
		Visualise3D.show(new Mapper() {

			@Override
			public double f(double x, double y) {
				// double h = 0.000001;
				// return (f.calculate(x, y) - f.calculate(x-h, y) )/h;
				// double[] val = f.calculateAll(x, y);
				// System.out.println(x+"\t"+y+"\t"+val[0]+"\t"+val[1]+"\t"+val[2]);
				// return val[1];//(x, y);//
				// double val = func.calculateAllFunctions(x,
				// y)[NodeList.indexOf(xi)].val;//*func.calculateAllFunctions(x,
				// y)[NodeList.indexOf(xi)+1].val;
				double val = func.calculate(x, y);
				// System.out.println(x+"\t"+y+"\t"+val+"\t"+f.calculate(x, y));
				return val;
				// return 0.0;
			}
		}, new Range(-0.5, 0.5), new Range(-0.5, 0.5));
	}

	public ArrayList<Double> getBaseVector(double x, double y, Particle xi) {

		return getBaseVector(x, y, xi.xi, xi.yi, xi.rho);
	}

	public ArrayList<Double> getBaseVector(double x, double y, double xi, double yi, double rho) {
		ArrayList<Double> L = new ArrayList<Double>();
		for (int i = 0; i <= this.PolinomialDegree; i++) {
			for (int j = 0; j <= this.PolinomialDegree; j++) {
				L.add(Math.pow((xi), i) * Math.pow(yi, j));
			}
		}
		// L.add(Math.exp(10*(xi.xi)));

		return L;
	}

	public ArrayList<Double> getBaseVectorDerivatyX(double x, double y, Particle xi) {
		return getBaseVectorDerivatyX(x, y, xi.xi, xi.yi, xi.rho);
	}

	public ArrayList<Double> getBaseVectorDerivatyX(double x, double y, double xi, double yi, double rho) {
		ArrayList<Double> L = new ArrayList<Double>();
		for (int i = 0; i <= this.PolinomialDegree; i++) {
			for (int j = 0; j <= this.PolinomialDegree; j++) {
				L.add(i * myPow((xi), i - 1) * myPow((yi), j));
			}
		}
		return L;
	}

	public ArrayList<Double> getBaseVectorDerivatyY(double x, double y, Particle xi) {
		return getBaseVectorDerivatyY(x, y, xi.xi, xi.yi, xi.rho);
	}

	public ArrayList<Double> getBaseVectorDerivatyY(double x, double y, double xi, double yi, double rho) {
		ArrayList<Double> L = new ArrayList<Double>();
		for (int i = 0; i <= this.PolinomialDegree; i++) {
			for (int j = 0; j <= this.PolinomialDegree; j++) {
				L.add(j * myPow((xi), i) * myPow((yi), j - 1));
			}
		}
		return L;
	}

	public List<Particle> getNodeList() {
		return NodeList;
	}

	public Particle getXi() {
		return xi;
	}

	public void setNodeList(List<Particle> nodeList) {
		NodeList = nodeList;
		// nodeList.remove(xi);
		// nodeList.add(0, xi);
		// for (Iterator<Particle> iterator = nodeList.iterator();
		// iterator.hasNext();) {
		// Particle particle = iterator.next();
		// if (particle==this.xi) {
		// nodeList.remove(particle);
		// nodeList.add(0, particle);
		// }
		//
		// }
	}

	private double myPow(double x, double pow) {
		if (pow < 0) {
			return 0.0;
		} else {
			return Math.pow(x, pow);
		}
	}

	protected void calculateAlfa(double x, double y, Particle xi) {
		ArrayList<Double> BaseVector = getBaseVector(x, y, xi);
		int n = BaseVector.size();
		for (int i = 0; i < n; i++) {
			double val = BaseVector.get(i) * func.calculate(x, y, xi.xi, xi.yi, xi.rho);
			for (int j = 0; j < n; j++) {
				A[i][j] += val * BaseVector.get(j);
			}
		}
	}

	public static void main(String[] args) {
		final ArrayList<Particle> nodes = new ArrayList<Particle>();
		// ????????? ???????:
		double ax = -0.5, bx = 0.5, ay = -0.5, by = 0.5;

		int n = 3, m = 3;

		double hx = (bx - ax) / (double) n, hy = (by - ay) / (double) m;

		// ?????+???????:
		for (int i = 0; i <= n; i++) {
			for (int j = 0; j <= m; j++) {
				Particle p = new Particle(ax + i * hx, ay + j * hy, 2.5 * hx);
				nodes.add(p);
			}
		}
		Gauss2D gau = new Gauss2D(20);

		// final BaseFunction ff = new BaseFunction(1,
		// nodes.get(nodes.size()/2), new SquareSupport(new CubicSpline()));
		final BaseFunction ff = new BaseFunction(2, nodes.get(0), new SquareSupport(new CubicSpline()));
		final BaseFunction ff1 = new BaseFunction(2, nodes.get(1), new SquareSupport(new CubicSpline()));
		// final BaseFunction ff = new BaseFunction(1, new Particle(0.5, 0.5,
		// 2.5*hx), new SquareSupport(new CubicSpline()));
		ff.setNodeList(nodes);
		ff1.setNodeList(nodes);
		ff.display();
		ff1.display();
		// 0.014842156360497365
		// 0.16056373301574658
		System.err.println(Math.abs(gau.integrate(-0.5, 0.5, -0.5, 0.5, new IFunction2D() {

			@Override
			public double calculate(double x, double y) {
				// double val = ff.calculate(x, y);
				// double val2 = ff1.calculate(x, y);
				Values[] v = ff.calculateAllFunctions(x, y);
				double val = v[nodes.indexOf(ff.xi)].val;
				double val2 = v[nodes.indexOf(ff1.xi)].val;
				return val * val2;
			}

			@Override
			public double calculateDerivatyX(double x, double y) {

				return 0;
			}

			@Override
			public double calculateDerivatyY(double x, double y) {

				return 0;
			}
		})));

		System.err.println(Math.abs(gau.integrate(-0.5, 0.5, -0.5, 0.5, new IFunction2D() {

			@Override
			public double calculate(double x, double y) {
				// double val = ff.calculateDerivatyX(x, y);
				// double val2 = ff1.calculateDerivatyX(x, y);
				long start = System.currentTimeMillis();
				Values[] v = ff.calculateAllFunctions(x, y);
				System.out.println(System.currentTimeMillis() - start);

				double val = v[nodes.indexOf(ff.xi)].valX;
				double val2 = v[nodes.indexOf(ff1.xi)].valX;
				return val * val2;
			}

			@Override
			public double calculateDerivatyX(double x, double y) {

				return 0;
			}

			@Override
			public double calculateDerivatyY(double x, double y) {

				return 0;
			}
		})));
	}

}
